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We present three-dimensional solutions of the magnetohydrostatic equations in the co-rotating 
frame of reference outside a magnetized rigidly rotating cylinder. We make no symmetry as- 
sumption for the magnetic field, but to be able to make analytical progress we neglect out- 
flows and specify a particular form for the current density. The magnetohydrostatic equations 
can then be reduced to a single linear partial differential equation for a pseudo-potential U, 
from which the magnetic field can be calculated by differentiation. The equation for U can 
be solved by standard methods. The solutions can also be used to determine the plasma 
pressure, density and temperature as functions of all three spatial coordinates. Despite the 
obvious limitations of this approach, it can for example be used as a simple tool to create 
three-dimensional models for the closed field line regions of rotating magnetospheres without 
rotational symmetry. 
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Introduction 



Three-dimensional analytical solutions of the magnetohydrostatic equations are 
difficult to find and only few solutions are known explicitly. A systematic method 
for finding a special class of three-dimensiona l magnetohydros t atic equilibria has 
been developed by in a s eries of papers by Low (|l985l . ll99ll . [T99l . ll993al . ll993bl l and 
Bogdan and Low ( 19861 ). The method relies on the presence of an external force 
derived from a potential, for example gravitation, and assumes that the electric 
current density has specific properties to allow analytical progress. In the simplest 
possible case, the fundamental equation to be solved is linear and it can be shown 
that in Cartesian and sp herical geometry with external gravity it is ve ry similar to 
a Schrodinger equation (jNeukirchl Il995l : iNeukirch and Rastatterl Il999l ) . Therefore, 
in the linear case, solutions can be found u sing standard methods, like expan- 



sion in terms of orthogonal function systems ( Rudenko 200 ll ) or Green's functions 
(IPetrie and Neukirch! j200ol ^ . Some solutions have also been found for nonlinear 



cases ( Neukirch! 1 1 99?! ) . Solutions found by this method have, for example , been 



used for models of solar and stellar corqnae (e .g. Zhao and Hoeksema 19931. 



Gibson and Bagenall[l995l : lGibson e7^ll99fil : IZhao et ai.ll2non[ 



Ruan et al. 



1994 : 



20081 : 



Lanzal 20081 ) ■ Another, but less general method for finding three-dimensional mag- 



netohysdrostatic equilibria in external gravitational fields has been proposed by 
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Osherovich (Il985al . Il985bl ) 



So far the method has only been used to find three-dimensional solutions of the 
magnetohydrosta tic equation s in Cartesian or spherical geometry with external 
gravity, although Lowl ( 1991 ) also generalizes the method to more general exter- 



nal potentials and forces, including in particular the centrifugal force in a rigidly 
rotating system. So far, no three-dimensional solutions to the magnetohydrostatic 
equations have been found for this case and it is the aim of this paper to close 
this gap. Such solutions could, for example, be used to model the closed field line 
regions of rotating magnetospheres without symmetries, such as those of the outer 
planets. 

We restrict our analysis to the somewhat artificial case of magnetic fields outside 
a rigidly rotating cylinder, but purely for reasons of mathematical convenience as 
it is much easier to deal with the boundary conditions on a cylindrical boundary. 
It is in principle also possible to treat the more realistic case of a spherical body, 
but imposing boundary conditions would be more difficult. 

Planetary magnetospheres are embedded in the solar wind (or a stellar wind 
in the case of exoplanets) and the boundary between the closed field line regions 
and the solar wind/open field line regions is a free boundary determined by total 
pressure balance, which would have to be calculated as part of the solution. The 
same statement applies to modeling the closed field line region of a stellar magne- 
tosphere, where the boundary between closed and open field line regions carrying a 
stellar wind would have to be determined as part of the solution (in a first idealized 
step one could maybe try to model the open field line regions by a potential field). 
There is little hope that any analytical progress could be made if the calculation 
of these free boundaries is included in the problem, as finding three-dimensional 
analytic solutions of the magnetohydrostatic equations is a formidable problem 
even without considering free boundaries. The solutions presented in this paper 
should therefore only be considered as a step towards solving the global problem 
of modeling magnetospheres, but not as the full solution. 

The paper is organised as follows. In section [2] we present the basic theoretical 
framework used in this paper, in section [3] we present a number of analytical three- 
dimensional solutions of the magnetohydrostatic equations in cylindrical geometry 
and we present a summary and conclusions in section HI 



2. Basic theory 



We consider a cylinder of radius R and infinite length rotating rigidly with constant 
angular velocity $7 about its symmetry axis, which we take to be identical with the 
z-axis. We will use a cylindrical coordinate system w, z in a frame of reference 

co-rotating with the cylinder. 

In this frame of reference the MHD equations are (see e.g. Mestel 19991 . chapter 

5) 



j X B - Vp - pSJV = 0, 
V X B = ;Uoj, 
V • B = 0, 



(la) 
(lb) 
(Ic) 



where B is the magnetic field, j is the current density, p is the pressure, p is the 
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plasma density and 

V = -^n^w^ (2) 

is the centri fugal poten tial. 
Following (|i991I ). we assume that 

^oj = VF X VV. (3) 

Substituting this into the force balance equation (jlap results in 

— (B • VF)VV - — (B • VV)VF -Vp- pVV = 0, (4) 



which implies that 



It follows that 



p(o7,0,i)=p(F,V). (5) 



''=-(^).^>-^^'' 

The possibi lity of furt her analytical progress now depends on the choice of the 
function F. ILo^ (|l99lh suggests the choice 

F{w, z) = k{V)^ ■ VI/, (7) 

where niV) is a free function of the theory. The current density is then given by 

m = k(f)v(b • vv) X vy, (8) 

and because the current density depends linearly on the magnetic field, Ampere's 
law (lib!) is also linear in B. Other choices for F are possible, but l ead to Anipere 's 
law being nonlinear in B (for an example in Cartesian geometry see NeukirchI 1997 ). 
Substituting ([7]) into (j6ap we can rewrite it in the form 

' F. (9) 



Direct integration gives 



P = Po(.V)-- —F\ (10) 



Substituting ([7]) for F into (llOp . we finally obtain 

1 



P = PoiV)-—K{V)iB-VVf (11) 
2/Uo 
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for the pressure. Here, po{V) is a positive function which represents a hydrostatic 
background pressure. 

The expression for the density can be evaluated by using (fTOl) to get 



dp\ _dpo 1 (Ik 2 .,2N 

dvjp dV 2fXoK^{V)dV ■ ^ ' 

Substituting this into (I6b[) and using (I7|) the final expression for the density finally 
becomes 

This form of the plasma pressure and density is a direct consequence of the assum- 
ing that current density has the form ([3]) with the function F given by ([7|). 
It can be seen that the term representing the background density, po{V), is given 

by 

MV) = (14) 

For the case discussed in the present paper, V is strictly negative and the density 
is strictly positive, which implies the background pressure must be monotonically 
increasing with distance from the rotation axis, with the precise form depending on 
the form of the background density. We shall not present any detailed discussion 
of a specific model for the hydrostatic background plasma in this paper, but it is 
important to choose the hydrostatic background plasma in such a way that it guar- 
antees that the full density and pressure are always positive, because that is not 
automatically guaranteed. Also the background pressure will become important 
if one wants to solve the full magnetospheric problem including the free bound- 
aries because continuity of total pressure across the free boundary is a necessary 
condition for an equilibrium solution. 

Although in the present paper we will only make use of the potential ([2]) we 
remark that in deriving the expressions for the pressure and the density, (jlip and 
(jl3p . no use has been made of the particular coordinate system or form of the 
potential V. It should be emphasized, however, that this does of course not imply 
that this theory can be generalized to differentially ro tating systems as that would 



violate Ferraro's iso-rotation theorem (jFerrard 119371 ). A permitted generalization 



of the potential ([2|) would be to add a gravitational potential, i.e. instead of having 
just a magnetized central body one would have a massive magnetized central body. 
The shape of the body and thus the form of the potential V is irrelevant for the 
validity of (jlip and (jl3p . as long as the gravitational potential is time-independent 
in the rotating frame of reference. 

An expression for the plasma temperature can be obtained if we assume that the 
plasma satisfies the equation of state of an ideal gas, 

(-) 

where R is the universal gas constant and /i is the mean molecular weight. 

The magnetic field can be calculated using Ampere's law ()lbp . The current den- 
sity can be written as 



^^oj = V X (FVV) 



(16) 
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SO that 

V X B = V X (FVV). (17) 

This equation can be integrated with the result 

B = VU + FVV, (18) 

where U is a free function. Using ^ to substitute in for F, we obtain 

B = VU + k{V){B ■ VV)VV. (19) 

Note that this is not yet the final expression for B as the right hand side still 
depends on B. Taking the scalar product of (jl9p with W and solving for B • V^, 
we obtain 



vu-vv 

1 - K{v){vvy 



B-Vy= , .,,,,,,^,,,2 - (20) 



The final expression for the magnetic field therefore is 



Again, expression (j2ip is completely independent of the form of V, bearing in mind 
the constraints discussed before. We also emphasize again that the specific form of 
(f2T|) is a direct consequence of ^ and ([7]). 
For the form ^ of the potential we find the components of B to be 

_ 1 9[/ 

^ ~ 1-k{V){V')^ 5^ ' ^^^^^ 

B^ = -% (22b) 
B. = - (22c) 

with 

From (j2ip and the equations for the components of the magnetic field, (j22ap . ()22bp 
and (j22cp . we see that C/ is a pseudo-potential, from which B can be determined 
by differentiation. We call it a pseudo-potential because the expression for B.^ is 
modified compared to that derived from a normal potential, where B would be 
given by the the gradient of the potential. This modification is due to the presence 
of currents in the system, represented by the factor 1/(1 — kV''^) in the equation 
for B^, (l22iD . 

Finally, the pseudo-potential U is determined by substituting (j2ip into the 
solenoidal condition (fTc|) . giving 



V . ( Vf/ + . VV)VV) = 0. (24) 
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Again, this equation is valid for general forms of the potential V , and not just for 
the form used in the present paper. An alternative form of this equation is 

V • (M • VU) = 0, (25) 

with the 3x3 matrix M defined as 



Here I is the 3x3 unit matrix. Equation (j25p is particularly useful if W has 
more than one non-vanishing component. This would, for example, be the case 
if we would treat the present problem in spherical polar coordinates instead of 
cylindrical coordinates. 
For the centrifugal potential ([2]) we define 



i{w) = ^'^w'^k{V), (27) 



and rewrite (IMl) as 



I d f w dU\ 1 d'^U d'^U , , 

+ -oTr7T + lT-T = 0. (28) 



w dw \l — ^{w) dw J w"^ dcj)^ dz 

Equation (|28|) is the fundamental equation that has to be solved to calculate the 
three-dimensional magnetic field. The functions k.{V), or alternatively ^(tu) are 
free functions and can be chosen to make analytical progress possible. We remark 
that it is only possible to choose directly the function S,{vj) instead of k{V) because 
of the one-to-one mapping between zu and V defined by ([2]). In cases where this 
mapping is not one-to-one, more care has to be taken. 

Using the definition (f27|) of .^(a;) the expressions for the pressure and the density 
can be rewritten as 

P = Po{w -— / — , 29a 

2/Uo (1 - 4(ro)) yomj 



1 



dpo I I ( duV 1 



dw 2/xo (1 — ^(ro))2 dro \dw J /xq 



--e(B-v)i?. 



(29b) 



An important point to note is that (j29a|) and ()29bp (or alternatively (fTT|) and 
(fTSl) ) do not automatically guarantee that pressure and density are always positive, 
but that one can always arrange them to be positive by choosing an appropriate 
background plasma model. 



3. Example Solutions 



Equation (I28p is very similar to Laplace's equation and one can immediately see 
that it has separable solutions of the form 



U{w, (j), z) = Fmk{'^) exp(im(^) exp(iA:z). 



(30) 
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with m integer and k real. If we substitute (I30p into (I28p we find that the radial 
function Fmk{'^) satisfies the equation 



w aw yi — E^yw) aw ) yw^ ) 

This ordinary second order differential equation will generally have two linearly in- 
dependent solutions, F^^^iw^ and F^l{w), say. Since the partial differential equa- 
tion for U is linear the solutions for different m and k may be superposed to give 
the general solution in the form 



U{w,(t>,z)= expiimcj)) dk exp{ikz)[Am{k)F^l{w) + Bm{k)F^l{w)]. 

m=-oo J-oo 

(32) 

The complex coefficients Am{k) and Bm{k) will be determined by the boundary 
conditions imposed on the surface of the cylinder {w = R) and on the outer 
boundary {w = Wout) and can be calculated using the standard formulae from the 
theory of Fourier transformations and series. From a mathematical point of view, 
the outer boundary could be at infinity, but from a physical point of view the 
problem will usually only be considered on a finite domain since the assumptions 
made for deriving our theory will break down at a finite distance from the rotation 
axis. 

For general choices of KiV)^ or alternatively, £,{w), (p8|) . or the corresponding 
radial equation (j3ip still has to be solved numerically. For some cases, however, 
exact solutions can be found, and we will discuss a couple of these cases in more 
detail now. In the following we assume that w and z have been normalised by 
the radius of the cylinder i?, i.e. that the cylinder radius is I in these normalised 
coordinates. 

The simplest possible choice for £,{w) is 

^{w) = ^0 = constant. (33) 
In this case Equation (I3ip takes the form 

1 d / w dFmk \ f , ,^2 



w aw \1 — qo aw J \ 



with the two linearly independent solutions 



F)nA^) = Iu{k^/l - Co ro), (35a) 
Fj^liw) = K,ik,/r^o rv), (35b) 



wher e I^ix) and Ky{x) are modified Bessel functions (jAbramowitz and Stegun 
196.51 ) with u = niy/l — Cq. Here < 1 has been assumed. 



Another function £,{w) for which analytical solutions of pip can be found is 



C{w) = 1 — qw'^, q = constant > 0. 



(36) 



55 



Geophysical and Astrophysical Fluid Dynamics tncentrircviscdedited 



8 T. Neukirch 

In this case (I3ip has the form 

,2 



^^(±^^)f^,,A,,„, = 0. (37) 

VJ aw \ qw azu J \ ■w^ J 



With the coordinate transformation 



w = ^tn^, (38) 



371) can be transformed into 



<5?Fmk (m? _ ,2 



\i^ + k'\ Fmk = 0. (39) 



The solutions of (j39p can be expr essed in terms of confluent hyp ergeometric func- 
tions M{a,b,x) and U{a,b,x) as ( Abramowitz and Stegun 19651 . Chapter 13), 



F^liw) = 2^/^zuexp{-^/^zu)M{l + m'^/{8k'^),l,2^/^w) (40a) 
F^liw) = 2^/^wexp{-^/^w)U{l + m'^/{8k^),l,2^/^w). (40b) 

To calculate model magnetospheres one would now have to choose appropriate 
boundary conditions, determine the coefficients Afn{k) and Bm{k) and calculate 
the pseudo-potential U using (j32p . From ?7, the magnetic field, the plasma pres- 
sure and plasma density can be derived using (j22ap to (|22cp and (j29ap to ()29bp . 
This will usually involve expressions containing an infinite series and an integral, 
which, although not inherently difficult, can make the algebra a bit tedious, espe- 
cially when calculating the pressure and density as terms that are quadratic in the 
magnetic field and its derivatives are involved. 

For the sake of simplicity we will use a different method to calculate and discuss 
an illustrative example solution for the case ^{w) = < 1- Defining 



28]) becomes 



^ = \/l - Co^^, (41a) 
4> = -J^<^^ (41b) 



\ d ( d\j\ 1 d'^u d^u ^ 

-TT^h^TT^ +— 7^ + 7rT = 0' 42 
w ow \ ow J ocp^ az"^ 

which is the Laplace equation in cylindrical coordinates w, (j) This shows 

that by rescaling of the coordinates and (p solutions to the Laplace equation 
H) can be transformed in to solutions of (j28p . A similar method has been used by 



Petrie and NeukirchI (|200d ) for finding a closed expression for a Green's function 



for the constant ^-case in Cartesian coordinates. The only constraint is that the 
transformed solutions have to be 27r-periodic in the (;/)-coordinate. It should be noted 
that this is not a requirement for the solutions of the transformed Laplace equa- 
tion which could be used to satisfy the periodicity condition in the untransformed 
coordinate (j). Another important point is that we are of course not restricted to 
using solutions to the Laplace equation found in cylindrical coordinates, instead 
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any other coordinate system could be used, e.g. spherical coordinates. We could 
even allow an inner boundary with a shape different from a cylinder, e.g. a sphere. 
The major disadvantage of this solution method compared to the standard sepa- 
ration of variables method is that it will be impossible to match given boundary 
conditions. Also, the method only works for the ^(nj) = constant case, but not for 
other choices of Ci'^)- 

Despite these limitations we shall now use the transformation method to gen- 
erate and discuss an illustrative solution of the general method described in this 
paper. As we are only interested in presenting an example solution we shall simply 
assume boundary conditions which match the solution we arrive at by using the 
transformation. As a starting point we shall use a magnetic dipole field but we 
allow the position of the dipole to be off-axis and the direction of the dipole to 
be non-aligned with the symmetry axis of the cylinder. Without loss of generality 
we can assume that the position of the dipole is = Wd, (p = and z = 0. In 
the corresponding Cartesian coordinate system, the dipole sits on the x-axis at a 
distance ro^ from the origin. 

The appropriate solution of the Laplace equation is 



- /io"i sin cos ^'(tz7 cos (/) — zurf) -|- sin sin ^'ro sin (/) -|- cos 0z 

U{W,(I),Z) = — — — ^ 7,-2, ' ^'^^^ 

47r (ro^ — 2ci7rfCi7 cos -|- roi + z^j'^/^ 



where m is the magnitude of the magnetic dipole moment, ^ is the angle between 
the direction of the dipole moment and the direction of the ro-axis and the angle 
between the direction of the dipole moment and the z-axis. 

In order to make the transformed solution 27r-periodic, the argument of the 
trigonometric functions in U has to be an integer multiple of (j). This can be satisfied 
by letting, for example, 



Co = I (44a) 
1 



\/Wo=2' (44b) 
^ = 24). (44c) 



Other possible choices for the value of are Co = 1 ~ n~^, n > 2 integer, which 
gives <f> = n(j). 

Normalising w and z by the cylinder radius R and the potential U by nom/ [ttB?), 
we obtain 



sin cos ^'(ro cos(2(^) — rod) -|- sin sin ^'Ci7 sin(2(/)) -|- 2 cos 02; 
U(w, </., z) = (^2 - 2vodVJ cos(2,/.) + wl+ Az^f/^ ' 

(45) 

where tu^ = 2wii- 
Using the abbreviations 



D = sin cos ^(tu cos(2(/)) — Wd) + sin sin '^■w sin(2(/)) -|- 2 cos 0z, (46a) 
r = (ro^ - 2wdV0 cos(2(/)) + wl + Az^ fl'^, (46b) 
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Figure 1. Three different viewing angles (<j> = 0, 4> = 2n/3 and <j> = Ait/3) of a 3D plot of selected field 
lines for the example solution for = = 3/4 and the dipole parameters zu^i = 0.5, = 5r/4 and 

* = 7r/4. 

the magnetic field components are given by 

l2D{w - WdCos{2(p)) - 4sin6(cos^cos(2(/») + sin ^ sin(2(/>))r2 



6Dwd sin(2(/)) - 2 sin 0(sin ^' cos(20) - cos sin(20))r^ 



UDz - 2 cos Qr^ 



, (47a) 
(47b) 
(47c) 



The normahsing factor for the magnetic field is Bq = {iiom) / {t: B?) . In Figured] 
we show the magnetic field for the case 'Wd = 0-5R, Q = 7r/4 and ^ = 7r/4. The 
same set of field lines is shown from three different perspectives corresponding to 
the viewing angles </> = 0, 2tt/3 and 47r/3. 

The pressure (normalised to Bq/{2^q) = {hqw?) / (2-k'^ R^)) is given by 



p{w,(l),z) =po{w) 

[3D{w - Wd cos(2(/))) - sin 0(cos cos(2(^) + sin ^' sm.{2(j)))r'^]'^ 



- 12^ 

and the density is given by 

p{w,(t),z) 



,10 



(48) 



1 
2w 



Aw 



(49) 



with the density normalised by i^Lom? /{-k^^R^). We have refrained from giving 
the explicit expression for the second term in the density as it is quite lengthy, 
but not really contributing to any better understanding of the problem. We are 
still free to specify a hydrostatic background plasma to complete our example. In 
Figure [2] we show contour plots in the 0-z-plane for w = 1.5 and w = 3.0 of 
the logarithm of the pressure and the density. The background terms have been 
subtracted to emphasize the variation of pressure and density in (p and z. It turns 
out that for the two radii shown, the pressure term is always negative and we 
therefore show contours of the natural logarithm of its modulus, whereas the density 
term is positive. This means that the plasma pressure will be reduced compared 
to whatever background pressure is chosen, whereas the plasma density will be 
enhanced above the background plasma density. One can also see the strong radial 
decrease of the density term by comparing the values of the contours at w = 1.5 
and w = 3. This is consistent with ()49p which shows that the density term should 
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og(lpl) \oq{\p I) 




Figure 2. Contour plots in the (p-z plane at radii = 1.5 (top row) and ro = 3.0 (bottom row) of the 
natural logarithm of the pressure and the density variation. Note that the background quantities have been 
subtracted as they can be chosen independently and only the three-dimensional terms are shown. 

have a Ci7~^-dependence, whereas the pressure term has a variation (see (j48p ). 
Although the fuh expressions for the density and pressure will depend on the 
background model, this will usually mean that the three-dimensional variations of 
pressure and density will only be of importance close to the central cylinder and 
decrease fast with distance from the rotation axis {w). 



4. Summary and Discussion 



Following [Lo3 (|l99lh we have presented a relatively simple analytical approach to 



calculate three-dimensional analytical MHD equilibria of rigidly rotating magne- 
tospheres in cylindrical geometry. The fundamental equation that has to be solved 
is a linear second order PDE for a pseudo-potential U. The theory contains a free 
function k{V), or alternatively a function ^(tn), which can be chosen such that 
analytical progress is possible. Using the standard method of separation of vari- 
ables, analytical solutions have been found for two different choices of (,{zu). For the 
choice of = constant, one can also find solutions using a coordinate trans- 
formation method. We have used this method to find and discuss some illustrative 
solutions. 

The solutions found here do not represent a complete model of the closed field 
line regions of rotating magnetospheres, because such a model would include the 
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calculation of the free boundaries between open and closed field line regions and/or 
the boundaries between the closed field line regions and the plasma in which the 
magnetosphere is embedded. Determining the free boundaries generally involves 
matching the total pressure on both sides of the free boundary and obviously the 
choice of the hydrostatic background plasma, which is part of the set of solutions 
presented here will be important. Generally, if one would attempt to calculate a 
model with free boundaries using the results presented in this paper the separation 
of variable method is better suited than the transformation method due to the 
availability of a complete set of functions allowing the expansion of the general 
solution. The difficulty would be the nonlinear nature of the boundary conditions 
on and the unknown shape of the free boundary. 

Possible modifications and extensions of the theory are to use spherical inner and 
outer boundaries or to use a combination of centrifugal and gravitational potential. 
This could, for exampl e, be used to mo del the closed field line regions of fast 
rotating stars (see e.g. Rvan et al\ 20051 ). In cylindrical geometry the modelling 



process will be very similar to that presented in the present paper, but with the 
additional complication that the combined potential 

V{w) = -]p}vj^ + 2GM ln(w/i?) (50) 

is not monotonic and thus cannot be mapped one-to-one to tu, because the com- 
bined potential has a maximum at the co-rotation radius where the Kepler 
velocity equals the rotational velocity VlWc- A simple combination of K;(y)y'^ into 
a function ^(lu) will usually lead to singularities in pressure and density at the 
co-rotation radius. Numerical solutions for physically sensible choices of k(V^) are, 
however, possible in both cylindrical and spherical geometry and will be the subject 
of a future publication. 
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